Coordinate - Based Meta - Analysis of fMRI Studies with R

This paper outlines how to conduct a simple meta-analysis of neuroimaging foci of activation in R. In particular, the first part of this paper reviews the nature of fMRI data, and presents a brief overview of the existing packages that can be used to analyze fMRI data in R. The second part illustrates how to handle fMRI data by showing how to visualize the results of different neuroimaging studies in a so-called orthographic view, where the spatial distribution of the foci of activation from different fMRI studies can be inspected visually. Functional MRI (fMRI) is one of the most important and powerful tools of neuroscientific research. Although not as commonly used for fMRI analysis as some specific applications such as SPM (Friston et al., 2006), AFNI (Cox and Hyde, 1997), or FSL (Smith et al., 2004), R does provide several packages that can be employed in neuroimaging research. These packages deal with a variety of topics, ranging from reading and manipulating fMRI datasets, to implementing sophisticated statistical models. The goal of this paper is to provide a brief introduction to fMRI analysis, and the various R packages that can be used to carry it out. As an example, it will show how to use simple R commands to read fMRI images and plot results from previous studies, which can then be visually compared. This is a special form of meta-analysis, and a common way to compare results from the existing literature. A brief introduction to fMRI data analysis Before getting into the details of the R code, it is important to have a clear picture of what type of data and analyses are common in fMRI research. An fMRI study consists of a set of 3D images that capture brain activity at fixed intervals of time while participants are performing a given task inside an MRI scanner. The interval between image acquisitions is known as repetition time or TR, and is typically 2 seconds. Each image is made of thousands of almost-cubic elements called voxels. The size of a voxel (typically on the order of 3× 3× 3 mm) represents the lower limit of the spatial resolution of an image. Not all the voxels are acquired at the same time: each image is typically created by acquiring each horizontal 2D plane (or “slice”) in serial order. Thus, voxels belonging to different slices are acquired at different times within the same 2-second TR. The MRI scanner captures brain “activity” only indirectly, by measuring the so-called Blood Oxygen-Level Dependent (BOLD) signal—that is, the amount of oxygenated blood present in every position in space. The BOLD signal is known to vary as a function of neural activity. Specifically, the BOLD response to an increase of neural activity at time t = 0 follows a specific time course, which is known as the hemodynamic response function h(t), and commonly modeled as the difference between two gamma functions (Worsley et al., 2002): h(t) = ( t d1 )a1 × e t−d1 b1 − c ( t d2 )a2 × e t−d2 b2 Figure 1 illustrates the typical shape of the hemodynamic response function, as returned by the function fmri.stimulus of the fmri package (Polzehl and Tabelow, 2007; Tabelow and Polzehl, 2011) using the following code. library(fmri) t <seq(0, 30, 0.1) # from 0 to 30 secs, in increments of 100ms h <fmri.stimulus(301, durations = c(0), # h(t) for instantaneous event at t=0, onsets = c(1), rt = 0.1), # Sampled at TR = 0.1 secs plot(t, h, type = "l", lwd = 2, col = "red", xlab = "Time (secs)", ylab = "h(t)", main = "Hemodynamic Response Function") Note how the function peaks only about 6 seconds after the event that triggered neural activity: The sluggishness of this response is one of the major drawbacks of functional neuroimaging, as it must be accounted for when designing an experiment and analyzing the data. After collection, fMRI images need to be preprocessed to remove various sources of noise. A typical preprocessing procedure (e.g., Friston et al., 2006) consists of the following steps: 1. Spatial realignment, whereby head motion from one image to another is removed by means of rigid-body transformation; The R Journal Vol. 6/2, December 2014 ISSN 2073-4859 CONTRIBUTED RESEARCH ARTICLES 6 Figure 1: The hemodynamic response function, as implemented in the fmri package. 2. Slice timing correction, whereby the temporal differences in the time at which voxels belonging to different 2D slices are removed; 3. Normalization, where individual brains are warped onto a common anatomical template by means of a non-linear transformation. This step is crucial to aggregate and compare data collected from anatomically different brains; 4. Spatial smoothing, which is typically performed with a 3D Gaussian kernel and is aimed at reducing noise and removing residual differences in individual anatomy between participants. Specialized R packages exists that perform the first and the last of these steps. Specifically, package RNiftyReg (Clayden, 2013) provides functions for realigning images in the NIFTI file format, while spatial smoothing is available in the AnalyzeFMRI package (Bordier et al., 2011). In addition, the fmri package implements an advanced smoothing procedure that adapts to the constraints of the underlying anatomical structures. To the best of my knowledge, and as confirmed by Tabelow and Polzehl (2011), there are currently no R packages that support the second and third preprocessing steps, which must still be performed using other software. Statistical analysis of fMRI data FMRI data is typically analyzed with a mass-univariate approach (Friston et al., 2006), where a single model is fitted to the timecourse of hemodynamic activity of each voxel independently. This is done by first creating a n × m matrix of regressors X, where the m columns represent the time courses of different experimental conditions and the n rows represent the discrete time points at which the images were acquired. To account for the sluggishness of the BOLD response (Figure 1), the matrix X is created by convolving each column of a matrix S describing onsets and durations of the experimental stimuli with the hemodynamic response function h, so that X = h ∗ S. Data from all the voxels Y of a single participant is then analyzed by fitting the linear model Y = Xβ+ e, and thus finding the values of β that minimize the residual sum of squares eTe for each voxel. However, the experimenter is typically interested in examining differences between conditions, rather than the estimates of the conditions per se. These differences are often referred to as “contrasts”, and can be expressed as the product between a contrast vector c (usually a binary vector), and the vector of estimates β, i.e., cTβ. Thus, the null hypothesis H0 that there is no difference between conditions corresponds to the case cTβ = 0. The result of a neuroimaging analysis is a set of three-dimensional images. Each image is made of tens of thousands of voxels, and each voxel contains the statistical value of a test of a contrast vector c. These parametric images are typically thresholded at a level that corresponds to a significant value of the corresponding statistical test, corrected for multiple comparisons (Friston et al., 2006). This procedure yields a set of 3D clusters of spatially adjacent voxels, whose activity is similarly affected by the experimental manipulations. While R does not provide many tools for preprocessing, there are many R packages that can be used to analyze preprocessed fMRI datasets. For instance, the fmri package contains functions to fit statistical models to fMRI datasets in the way described above. In particular, the design matrix X is generated using the function fmri.design, and a linear model can be fit using fmri.lm. In The R Journal Vol. 6/2, December 2014 ISSN 2073-4859 CONTRIBUTED RESEARCH ARTICLES 7 addition, both the AnalyzeFMRI and the fmri packages also contain functions to perform Independent Component Analysis in both the spatial and temporal domains, as well as several procedures for multiple-error correction, such as Bonferroni, False Discovery Rate, and various instances of familywise corrections. The package RfmriVC (Bothmann and Kalus, 2013) permits users to modulate the effects of the hemodynamic response function with other variables that vary within a single session, such as time or fatigue. The package arf3DS4 (Weeda et al., 2011) uses the estimated βor t-images to identify active brain regions and their patterns of functional connectivity. The package BHMSMAfMRI (Sanyal and Ferreira, 2014) implements Bayesian multi-subject analysis. Finally, the package neuRosim (Welvaert et al., 2011) contains functions that simulate brain imaging data for testing purposes. A comprehensive list of these and related packages is maintained in the CRAN task view on Medical Image Analysis (http://CRAN.R-project.org/view=MedicalImaging). Describing the functionalities of these packages is beyond the scope of this paper; the interested reader can check the recent and extensive review by Eloyan et al. (2014) on the subject. Instead, the remainder of this article will introduce some simple code lines that demonstrate how to read imaging data, plot them over a brain image, and visually compare results from different neuroimaging studies. Incidentally, this type of comparison is a common way to familiarize oneself with a new field in neuroimaging research, and is quite often an important step in justifying hypotheses about the functions of one brain region. An example of coordinate-based meta-analysis To understand the coordinate-based meta-analysis, we need to take a step back. Remember that the result of an fMRI analysis is a set of clusters of voxels that are jointly above a certain statistical threshold. Because a single test might identify multiple clusters, neuroimaging papers typically focus only on those that are deemed the most interesting, and relegate the full list of clusters to separate tables. These tables detail various characteristics of each cluster, including its size (i.e., the number of voxels it spans) and position within the brain. To facilitate comparisons across studies, cluster positions are given in terms of three-dimensional coordinates in predefined stereotactic spaces (which will be discussed in the next section). Because a cluster spans many voxels, different conventions can be applied to determine which coordinates best represent a cluster’s position. For instance, some authors might report the coordinates of the cluster’s “peak” (i.e., the highest parameter value), while others report the cluster’s geometrical center or its center of mass. Fortunately, these values tend not to differ much from each other, and for simplicity we can consider them as functionally equivalent. In theory, and assuming they are available and comparable, one could directly enter contrast or statistical parameter images from different studies into a meta-analysis. This procedure is known as image-based meta-analysis, and is generally considered the most reliable (Salimi-Khorshidi et al., 2009). However, because cluster sizes are dependent on the statistical error correction procedures that were adopted, most meta-analyses of fMRI data, like the Activation Likelihood Estimate (Turkeltaub et al., 2002), focus instead on the consistency of cluster positions, i.e., how often a given region is reported, and how the spatial distribution of cluster peaks is affected by the study parameters. This procedure is known as coordinate-based meta-analysis, and has been successfully used to identify possible subdivisions within the same anatomical brain regions (for example, the rostral prefrontal cortex, Burgess et al. 2007; or the anterior cingulate cortex, Bush et al. 2000), or to suggest the involvement, in a certain cognitive function, of a brain structure that had been previously overlooked (e.g., the basal ganglia in language switching, Stocco et al. 2014). The remainder of this paper will introduce the packages and the R code that will be used to visualize foci of activation on a standard brain template. The final result will be an image like the one in Figure 2. Although simple, the code provided herein is useful in its own right. Images like Figure 2 have been published in many scientific papers, such as Bush et al. (2000) and Burgess et al. (2007). This very same code, in fact, has been used to generate two figures published in Stocco et al. (2014) and Buchweitz and Prat (2013). Displaying the MNI template The spatial coordinates of a cluster of voxels are given in reference to an ideal brain template, oriented in a given stereotactic space. Two stereotactic spaces are commonly used, the Talairach-Tournoux and the Montreal Neurological Institute (MNI) spaces. In this paper, we will assume that all the coordinates are in the MNI space, which is preferred by the majority of researchers (Poldrack et al., 2011). Talairach-Tournoux coordinates can be converted to their MNI equivalent using non-linear transformations, such as those proposed by Brett et al. (2001). The R Journal Vol. 6/2, December 2014 ISSN 2073-4859 CONTRIBUTED RESEARCH ARTICLES 8 Figure 2: A visual summary of the results of four neuroimaging studies. In the MNI space the x, y, and z axes correspond to the left-right, front-back, and up-down orientation of a brain. The origin of the space is located in a point in the middle of the brain, so that the brain volume spans both positive and negative coordinates along all axes. The first step to producing Figure 2 is to visualize the MNI template brain on an R plot. Many versions of the MNI template exist; in this paper, I am going to use the so-called "Colin27" version (Holmes et al., 1998). This template is available on various websites, and is included in many common software packages for fMRI, such as MRICron (Rorden et al., 2000). This example will use the template that is available at http://packages.bic.mni.mcgill.ca/mni-models/colin27/mni_colin27_1998_ nifti.zip. The template is encoded in a particular file format known as NIfTI, which is, by far, the most common file format for functional neuroimages. All the packages mentioned above contain functions to read and write NIfTI files. In this example, I will use the oro.nifti package (Whitcher et al., 2011), which has been developed specifically for this purpose. Assuming that the template’s URL is stored in the variable colin.url, the image can be easily loaded with this code: library(oro.nifti) temp <file.path(tempdir(), "mni_colin27_1998_nifti.zip") colin.url <"http://packages.bic.mni.mcgill.ca/mni-models/colin27/mni_colin27_1998_nifti.zip" download.file(colin.url, dest = temp) unzip(temp, exdir = tempdir()) colin <readNIfTI(file.path(tempdir(), "colin27_t1_tal_lin.nii")) The colin object is simply a 3D matrix. The matrix dimensions are of 181× 217× 181 voxels, and can be accessed by calling the dim function: dim(colin). Each slice of this matrix is in itself a 2D image, and can be visualized with any standard 2D plotting tools in R. For instance, a professional-looking sagittal (i.e., lateral) view of the brain can be visualized by calling the image function on the 80th 2D slice across the x axis, e.g., image(colin[80, , ], col = grey(0:255 / 255)), as shown in Figure 3. Notice that the original template, as shown in Figure 3, contains the brain as well as a surrounding skull. For clarity purposes, the skull has been omitted in Figure 2. The Colin27 archive that had been downloaded in the previous step also contains a “brain mask”, i.e., a binary image that spans the entire brain but leaves out the skull. An skull-stripped image of the brain can be generated as the point-wise product of the colin image and the mask: mask <readNIfTI(file.path(tempdir(), "colin27_t1_tal_lin_mask.nii")) colin <colin * mask The R Journal Vol. 6/2, December 2014 ISSN 2073-4859 CONTRIBUTED RESEARCH ARTICLES 9 Figure 3: The 80th 2D sagittal slice of the colin template. The three axes of the colin object are already aligned with the three axes of the colin template. However, the coordinate system of the image and the coordinate system of the colin image have different centers. In the colin object, the point of coordinates x = 0, y = 0, z = 0 is located at the leftmost voxel of the first row of the bottom slice. In the MNI space, however, the origin is located in the middle of the brain. The MNI origin corresponds to the voxel of coordinates x = 91, y = 126, z = 72 in the image space. Thus, the MNI coordinates need to be translated in space to be visualized correctly. Because each voxel in colin is a 1× 1× 1 mm cube, this affine translation can be easily implemented using sweep: # Origin of MNI space origin <c(x = 91, y = 126, z = 72) # Converting back and forth to MNI space mni2xyz <function(x) sweep(x, 2, as.double(origin), "+") xyz2mni <function(x) sweep(x, 2, as.double(origin), "-") Figure 2 is a rendition of the brain in so-called orthographic view, which is particularly helpful for visualizing the relative position of foci within the brain. The orthographic view consists of three orthogonal brain slices that intersect in a specific point. We will refer to this point as the center C of the orthographic view, with coordinates xC, yC, zC. The three views are produced by simply selecting the 2D slice of the 3D image with the corresponding values of C. The top left view in Figure 2 is called coronal view, and corresponds to the case where y = yC. The top right view in Figure 2 is called sagittal view, and corresponds to the case where x = xC. Finally, the bottom-left view is known as axial view, and corresponds to the case where z = zC. The package oro.nifti has its own method (orthographic) to visualize a brain in orthographic view. In this example, however, we will implement the procedure from scratch, as it is fairly easy and allows for more control on the proper visualization. The first step towards producing the image in Figure 2 is to create a the plot window and divide it into four quadrants. R’s default layout method is a convenient way to do so: layout(matrix(1:4, ncol = 2, byrow = TRUE), heights = c(181, 217), widths = c(181, 217)) layout.show(4) The heights and widths parameters in the code are used to assign non-uniform heights and widths to the four quadrants. This is needed, in turn, because the colin template is not cubic, and it spans 181 voxels along the x and z axis, but 217 on the z axis. Next, we need to define our center point C. Initially, I will simply assume that the center of the orthographic view is also the origin in MNI space, i.e., the point of coordinates x = 0, y = 0, z = 0. The coronal, sagittal, and axial views of the orthographic figure can be simply generated by subsequent The R Journal Vol. 6/2, December 2014 ISSN 2073-4859 CONTRIBUTED RESEARCH ARTICLES 10 calls of the default method image: center <origin greys <grey(1:max(colin) / max(colin)) par(mar = c(2, 2, 2, 2)) image(colin[, center["x"], ], col = greys) # Coronal image(colin[center["y"], , ], col = greys) # Sagittal image(colin[, , center["z"]], col = greys) # Axial plot.new() # Empty quadrant In the code above, the greys variable is used to create as many shades of grey as there are values of intensity in the Colin27 template. These colors as passed as an argument to visualize the brain slices in the canonical greyscale (instead of R’s default heat colors) like in in Figure 3. The mar variable sets the plot margins for each image to 2 lines of text, which is the minimum amount necessary to visualize the image axes in each view. Finally, the call to plot.new is required because R expects one plot for each quadrant defined by layout.